library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
library(opentripplanner)
library(tidyverse)
## -- Attaching packages ------------------------------------------------------------------------- tidyverse 1.3.0 --
## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.2
## v tidyr   1.1.2     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0
## -- Conflicts ---------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(ggthemes)
library(ggspatial)
library(smoothr)
## 
## Attaching package: 'smoothr'
## The following object is masked from 'package:stats':
## 
##     smooth
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Loading data and cleaning them

MA_state_plane <- "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs"

WGS84 <- "+proj=longlat +datum=WGS84 +no_defs"

subway_stops <- st_read("https://maps-massgis.opendata.arcgis.com/datasets/a9e4d01cbfae407fbf5afe67c5382fde_0.kml")

neighborhood_boundaries <- st_read("http://bostonopendata-boston.opendata.arcgis.com/datasets/3525b0ee6e6b427f9aab5d0a1d0a1a28_0.kml?outSR=%7B%22latestWkid%22%3A2249%2C%22wkid%22%3A102686%7D")

northend_neighborhood <- neighborhood_boundaries[neighborhood_boundaries$Name %in% c("North End", "Downtown", "West End","Beacon Hill", "Chinatown", "Leather District"),] %>%
  st_transform(crs = MA_state_plane)

subway_stops <- subway_stops %>%
  st_transform(crs = MA_state_plane)

neighborhood_boundaries <- neighborhood_boundaries %>%
  st_transform(crs = MA_state_plane) %>%
  st_buffer(dist = 0)


boston_stops <- subway_stops[neighborhood_boundaries, ] %>%
  st_transform(crs = MA_state_plane)

northend_stops <- boston_stops[northend_neighborhood,] %>%
  st_transform(crs = MA_state_plane)

northend_neighborhood <- northend_neighborhood %>%
  st_transform(crs = MA_state_plane)

Data

The data I will be using is the neighborhood map of Boston and all the MBTA subway stops within Boston. I will also look at the Downtown area more closely (includes the neighborhoods: North End, West End, Downtown, Beacon HIll, Chinatown and the Leather District.

plot1 <- ggplot(neighborhood_boundaries)+
  geom_sf()+
  geom_sf(data = boston_stops)+
  theme_map()+
  labs(title = "Boston")

plot2 <- ggplot(northend_neighborhood)+
  geom_sf() +
  geom_sf(data = northend_stops, aes(fill = "Metro Stops"))+
  guides(fill = "legend")+
  theme_map()+
  theme(legend.position = c(-0.3,0.2), legend.text = element_text(size = 10))+
  labs(title = "Downtown Boston area", fill = "")

grid.arrange(plot1, plot2, nrow = 1)

Getting the OpenStreetMaps street data for Boston

boston_streets <- opq(bbox = 'Boston MA USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_sf()

boston_streets <- boston_streets$osm_lines %>%
  st_transform(crs = MA_state_plane)

ggplot(boston_streets) +
  geom_sf() +
  theme_map()

Connecting to Open Trip Planner

# path_otp <- otp_dl_jar("OTP")

path_data <- file.path(getwd(), "OTP")
path_otp <- paste(path_data, "otp.jar",sep = "/")

otp_build_graph(otp = path_otp, dir = path_data, memory = 1024)

otp_setup(otp = path_otp, dir = path_data, memory =1024)
otpcon <- otp_connect()

Getting data on 10-minute walk isochrones from all the MBTA subway stops in Boston

boston_stops <- boston_stops %>%
  st_transform(crs = WGS84)

iso_10min_walk_boston <-
  otp_isochrone(otpcon = otpcon, fromPlace = boston_stops,
                mode = "WALK", cutoffSec = 600) %>%
  st_transform(crs = MA_state_plane)

Getting data on 10-minute walk isochrones from the MBTA subway stops in the Downtown area

northend_stops <- northend_stops %>%
  st_transform(crs = WGS84)

iso_10min_walk_northend <-
  otp_isochrone(otpcon = otpcon, fromPlace = northend_stops,
                mode = "WALK", cutoffSec = 600) %>%
  st_transform(crs = MA_state_plane)

Getting data on 10-minute bike isochrones from the MBTA subway stops in the Downtown area

iso_10min_bike_northend <-
  otp_isochrone(otpcon = otpcon, fromPlace = northend_stops,
                mode = "BICYCLE", cutoffSec = 600) %>%
  st_transform(crs = MA_state_plane)

otp_stop()

Isochrone 1-A

This map shows all the 10-min walk isochrones from all the MBTA subway stops in Boston.

northend_stops <- northend_stops %>%
  st_transform(crs = MA_state_plane)

boston_stops <- boston_stops %>%
  st_transform(crs = MA_state_plane)

ggplot()+
  annotation_map_tile(type = "cartolight", zoomin = 0, progress = "none") +
  geom_sf(data = neighborhood_boundaries, alpha = 0.2)+
  geom_sf(data = iso_10min_walk_boston, fill = "#F0B791", alpha = 0.5)+
  theme_map()+
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

Isochrone 1-B

This map shows all the 10-min walk isochrones as one aggregated area within the Boston boundaries.

iso_boston_union <- st_intersection(neighborhood_boundaries, iso_10min_walk_boston) %>%
  st_union()

ggplot()+
  annotation_map_tile(type = "cartolight", zoomin = 0, progress = "none") +
  geom_sf(data = neighborhood_boundaries, alpha = 0.2)+
  geom_sf(data = iso_boston_union, fill = "#F0B791", alpha = 0.5)+
  theme_map()+
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

Isochrone 2-A

This map shows all the 10-min walk isochrones from the MBTA subway stops in the Downtown area.

northend_stops <- northend_stops %>%
  st_transform(MA_state_plane)

ggplot()+
  annotation_map_tile(type = "cartolight", zoomin = 0, progress = "none") +
  geom_sf(data = northend_neighborhood, alpha = 0.25)+
  geom_sf(data = iso_10min_walk_northend, fill = "#F0B791", alpha = 0.3)+
  geom_sf(data = northend_stops) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

Isochrone 2-B

This map shows the 10-min walk isochrones as one aggregated area within the Boston boundaries.

iso_northend_union <- st_intersection(northend_neighborhood, iso_10min_walk_northend) %>%
  st_union

ggplot()+
  annotation_map_tile(type = "cartolight", zoomin = 0, progress = "none") +
  geom_sf(data = northend_neighborhood, alpha = 0.25)+
  geom_sf(data = iso_northend_union, fill = "#F0B791", alpha = 0.5)+
  geom_sf(data = northend_stops) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

Isochrone 3-A

This map shows the 10-minute bike isochrones from the MBTA subway stops in the Downtown area. Not super helpful since there are a lot of stops in a small area and they all overlap.

ggplot()+
  annotation_map_tile(type = "cartolight", zoomin = 0, progress = "none") +
  geom_sf(data = northend_neighborhood, alpha = 0.25)+
  geom_sf(data = iso_10min_bike_northend, fill = "#7EA285", alpha = 0.1)+
  geom_sf(data = northend_stops) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

Isochrone 2-B

This map shows the 10-min bike isochrones as one aggregated area within the Boston boundaries. The entire area is filled in, showing that every part of the Downtown area is within a 10 minute bike ride of a subway stop.

iso_northend_bike_union <- st_intersection(northend_neighborhood, iso_10min_bike_northend) %>%
  st_union

ggplot()+
  annotation_map_tile(type = "cartolight", zoomin = 0, progress = "none") +
  geom_sf(data = northend_neighborhood, alpha = 0.25)+
  geom_sf(data = iso_northend_bike_union, fill = "#7EA285", alpha = 0.3)+
  geom_sf(data = northend_stops) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

Graph #1

This scatterplot shows the Downtown area neighborhoods and plots them along how many MBTA subway stops are located within each neighborhood and how much of the neighborhood is covered by the 10-minute walk isochrones. They are also distinguished in size by the size of the neighborhoods.

northend_neighborhood$area <- northend_neighborhood %>%
  st_area() %>%
  as.numeric()

northend_neighborhood <- northend_neighborhood %>%
  mutate(iso_area = NA)

for (i in c(1:nrow(northend_neighborhood))){
  northend_neighborhood$iso_area[i] = iso_northend_union %>%
    st_intersection(northend_neighborhood[i,]) %>%
    st_area()
}

northend_neighborhood <- northend_neighborhood %>%
  mutate(iso_percent = 
           northend_neighborhood$iso_area/northend_neighborhood$area * 100,
         num_stops = lengths(st_covers(geometry, northend_stops)))

ggplot(northend_neighborhood, aes(x = num_stops, y = iso_percent))+
  geom_point(aes(size = area)) +
  geom_text(aes(label = Name), hjust = 0.5, vjust =2, size = 3)+
  scale_size_continuous(name = "Area of\nneighborhoods (m^2)") +
  labs(x = "Number of subway stops", 
       y = "Percent of neighborhood area\ncovered by isochrones")+
  theme_minimal()